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ABSTRACT 



o 
o 

■ We present results from a series of cosmological SPH (smoothed particle hydrody- 
namics) simulations coupled with the P^M (Particle-Particle-Particle-Mesh) solver for 
the gravitational force. The simulations are designed to predict the statistical proper- 
ties of X-ray clusters of galaxies as well as to study the formation of galaxies. We have 

I seven simulation runs with different assumptions on the thermal state of the intracluster 

gas. Following the recent work by Pearce et al., we modify our SPH algorithm so as 
to phenomenologically incorporate the galaxy formation by decoupling the cooled gas 
^ ' particles from the hot gas particles. All the simulations employ 128^ particles both for 

■ dark matter and for gas components, and thus constitute the largest systematic cata- 
logues of simulated clusters in the SPH method performed so far. These enable us to 
compare the analytical predictions on statistical properties of X-ray clusters against our 
direct simulation results in an unbiased manner. We find that the luminosities of the 
simulated clusters are quite sensitive to the thermal history and also to the numerical 
resolution of the simulations, and thus are not reliable. On the other hand, the mass - 
temperature relation for the simulated clusters is fairly insensitive to the assumptions 
of the thermal state of the intracluster gas, robust against the numerical resolution, and 
in fact agrees well with the analytic prediction. Therefore the prediction for the X-ray 
temperature function of clusters on the basis of the Press - Schechter mass function 
and the virial equilibrium is fairly reliable. 



Subject headings: galaxies: clusters: general ~ cosmology: miscellaneous ~ methods: 
numerical 
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1. Introduction 

Clusters of galaxies have been widely used as probes to extract the cosmological information. 
In particular, the cluster abundances including the X-ray temperature function (XTF), X-ray lu- 
minosity function (XLF) and number counts turn out to put strong constraints on the cosmological 
density parameter f^o and the fluctuation amplitude dg (Henry & Arnaud 1991; White, Efstathiou 
& Frcnk 1993; Jing & Fang 1994; Viana & Liddle 1996; Eke, Cole & Frenk 1996; Kitayama & Suto 
1996, 1997; Kitayama, Sasaki & Suto 1998). Since the theoretical predictions for those purposes 
are usually based on the Press - Schechter mass function and the simple virial equilibrium model 
for the X-ray clusters, the reliability of the resulted constraints is heavily dependent on the validity 
of those assumptions in the observed X-ray clusters. In fact there are numerous realistic physical 
processes which would somehow invalidate the simplifying assumptions employed in the above pro- 
cedure; one-to-one correspondence between a virialized halo and an X-ray cluster, non-sphericity, 
substructure and merging, galaxy formation, radiative cooling, heating by the UV background and 
the supernova energy injection, etc. 

While these can be addressed by hydro dynamical simulations in principle, the relevant simu- 
lations are quite demanding. This is why most previous studies simply check the Press - Schechter 
mass function against the purely A'^-body simulations in cosmological volume (Suto 1993; Ueda, 
Itoh k, Suto 1993; Lacey & Cole 1994), or focus on hydrodynamical simulations of individual clus- 
ters in a constrained volume (Eke, Navarro & Frcnk 1998; Suginohara & Ostriker 1998; Yoshikawa, 
Itoh & Suto 1998); most previous hydrodynamical simulations of clusters in cosmological volume, 
on the other hand, lacked the numerical resolution to address the above question, and/or ignored 
the effect of radiative cooling (Evrard 1990; Cen & Ostriker 1992, 1994; Kang et al. 1994; Bryan & 
Norman 1998; Frenk et al. 1999) except for a few latest state-of-the-art simulations (Pearce et al. 
1999a,b; Cen & Ostriker 1999). 

In this paper, we present results from a series of cosmological SPH (smoothed particle hydro- 
dynamics) simulations in cosmological volume with sufficiently high resolution. They enable us to 
compute the statistical properties of simulated clusters in an unbiased manner, and thus to directly 
address the validity of the widely used model predictions for the mass — temperature relation and 
the resulting XTFs. 



2. Simulations 

Our simulation code is based on the P^M-SPH algorithm, and the P^M part of the code 
has been extensively used in our previous N-hody simulations (Jing k, Fang 1994; Jing & Suto 
1998). All the present runs employ Adm = 128^ dark matter particles and the same number of 

gas particles. We use the spline (S2) softening function for gravitational softening (Hockney &: 

1/3 

Eastwood 1981), and the softening length egrav is set to L\)oy_/{lQN-^y^), where -Lbox is the comoving 
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size of the simulation box. We set the minimum of SPH smoothing length as hmm = egrav/4, and 
adopt the ideal gas equation of state with an adiabatic index of 5/3. 

We consider a spatially-flat low-density CDM (cold dark matter) universe with the mean mass 
density parameter Qq = 0.3, cosmological constant Aq = 0.7, and the Hubble constant in units of 
100 km s^^ Mpc~^, h = 0.7. The power-law index of the primordial density fluctuation is set to 
n = 1. We adopt the baryon density parameter Jib = 0.015/i~^ (Copi et al. 1995) and the rms 
density fluctuations on Sh"^ Mpc scale, ag = 1.0, according to the constraints from the COBE 
normalization and cluster abundance (Bunn & White 1997; Kitayama & Suto 1997). 

The initial conditions for particle positions and velocities at z = 25 are generated using the 
COSMICS package (Bertschinger 1995). We prepare two different initial conditions for Lbox = 
150/i~^Mpc and 75/t~^Mpc to examine the effects of the numerical resolution. For the two different 
boxsizes, we consider three cases for the ICM (intracluster medium) thermal evolution (Table 1); 
L150A and L75A assume non-radiative evolution of the ICM (hereafter non-radiative models). 
L150C and L75C include the efi:ect of radiative cooling (pure cooling models). In L150UJ and 
L75UJ (multiphase models), both radiative cooling and UV-background heating are taken into 
account, and an additional modification of the SPH algorithm is implemented in order to avoid 
the artiflcial overcooling, following Pearce et al. (1999a,b). In addition, we have L150CJ which 
is identical to L150UJ but without the UV background. While the L75C model is evolved until 
z = 0.5, the simulations of the remaining six models are computed up to z = 0. 

The radiative cooling and the UV-background heating take account of the photoionization of 
Hi, He I and Hcll, the collisional ionization of Hi, Hel and Hell, the recombination of Hll, Hell and 
Hem, the dielectronic recombination of Heii, Compton cooling and the thermal bremsstrahlung 
emission. The UV-background flux density J{i') is assumed to be parameterized as 

J{i^) = J2i{z) (^) X 10-2^ [erg/s/cmVsr/Hz], (1) 

where i/^ is the Lyman-limit frequency, and we adopt the redshift evolution: 

= (2) 

following Vedcl et al. (1994). We use the semi-implicit scheme in integrating the thermal energy 
equation (Katz, Weinberg &; Hernquist (1996)). 

It is known that SPH simulations with the effect of radiative cooling produce unacceptably 
dense cooled clumps. While Suginohara & Ostriker (1998) ascribed this to some missing physical 
processes such as heat conduction and heating from supernova explosions, Thacker et al. (1998) and 
Pearce et al. (1999a) showed that this overcooling could be due to an artificial overestimate of hot 
gas density convolved with the nearby cold dense gas particles, and can be largely suppressed by 
simply decoupling the cold gas (T < lO'^K) from the hot component. This numerical treatment can 
be interpreted as a phenomenological prescription for the galaxy formation. Therefore following 
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their spirit, we adopt a slightly different approach — using the Jeans criterion for the cooled gas: 

where h is the smoothing length of gas particles, Pgas the gas density, Cg the sound speed, and G 
the gravitational constant. In practice, however, we made sure that the above condition is almost 
identical to the condition of T < lO^K adopted by Pearce et al. (1999a). Nevertheless we prefer 
using the Jeans criterion because it is rephrased in terms of the physical process. Apart from the 
fact that such cooled gas particles are ignored when computing the SPH density of hot gas particles, 
all the other SPH interactions are left unchanged. 



3. Results 

3.1. Cluster Identification and Radial Profiles 

At each epoch, we identify gravitationally bound objects using an adaptive "friend-of-friend" 
algorithm (Suto et al. 1992), and select objects with more than 200 dark matter and gas particles 
as clusters of galaxies. Specifically we use the SPH gas density in assigning the local linking length 
with the linking parameter of 0.5. With this parameter, the mass function of the selected clusters 
is nearly identical to the conventional "fricnd-of-friend" with the fixed linking length of 0.2 times 
the mean particle separation. The virial mass Al for each cluster is defined as the total mass at 
the virial radius within which the averaged mass density is ~ 187:'^Q.q'^Pc{z) predicted from non- 
linear spherical collapse model, where pdz) is the critical density of the universe at z (White, 
Efstathiou &; Prenk 1993; Kitayama &; Suto 1996). The X-ray luminosity is computed on the basis 
of the bolometric and band limited thermal bremsstrahlung emissivity (Rybicki & Lightman 1979), 
which ignores metal line emission. We also compute the mass weighted and emission weighted 
temperatures, T™ and Tx, using 2-10 keV band emission. For models L150UJ and L75UJ, we do 
not include contribution of those cold particles which satisfy the criterion (3) in computing Tx, 

and L, since they are not supposed to remain as a gaseous component if an appropriate star 
formation scheme is implemented in our simulation. 

In Figure 1, we show the spherically averaged profiles of dark matter and gas densities, ICM 
temperature, and cumulative X-ray luminosity for the most massive cluster in models L150A, 
L150UJ and L150C. Those for models L75A, L75UJ and L75C arc also depicted in Figure 2. The 
centers of clusters arc defined as the peaks of gas density. The profiles of dark matter are fairly 
well-approximated by the model proposed by (Navarro, Frcnk, & White 1997). Due to an artificial 
cooling catastrophe, the pure-cooling model (L150C and L75C) produces unacceptably high X-ray 
luminosity, consistent with Suginohara &; Ostriker (1998). This artificial feature can be largely 
removed by decoupling the cold gas from hot component, and models L150UJ and L75UJ yield a 
fairly reasonable range of X-ray luminosities. Nevertheless we do not claim that the luminosities 
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in L150UJ and L75UJ are reliable; the current modified scheme is to be interpreted as a tentative 
and phenomenological remedy at best and should be replaced by more realistic scheme for galaxy 
formation. Lewis et al. (1999) have obtained a very similar result based on a more sophisticated 
treatment of galaxy formation in their simulation. They show that when both cooling and star 
formation are included in the simulation, cluster luminosities increase only moderately (a factor 
3) compared to the non-radiative case; this amount of increase is significantly less than that in 
the pure-cooling case. In addition, wc find that the values of the luminosities are still significantly 
affected by the numerical resolution (§3.3). 

In Figures 3 and 4, we show the spherically averaged profiles of local dynamical timescale 
t^yn, 2-body heating timescale t2body (Steinmetz & White 1997) and cooling timescale tcooi for 
representative clusters with 3 different mass scale for models L150UJ and L75UJ, respectively. 
Each timescale is defined as 

^ 7^' ^"^^ 

t2body = V^G^MoMpLlnA' 

, 377.gasA:BT' , , 

- 2Aeooi(n,r)' 

where G is the gravitational constant, pt^t the total mass density , (Jid the 1-dimensional velocity 
dispersion of dark matter particles, In A the Coulomb logarithm, M^m the mass of a dark matter 
particle, Pdm the mass density of dark matter, n^^^ the number density of gas, the Boltzmann 
constant and Acooi{n,T) the cooling function. The expression for i2body is given in Steinmetz &: 
White (1997) and we set the value of the Coulomb logarithm to InA = 5 as a nominal value. 

For the clusters with M > 2 x 10^^ Mq in the models with Lbox = 150/i~^Mpc , tcooi is shorter 
than t2body) and hence the artificial 2-body heating is not effective for these clusters. On the other 
hand, clusters with M < IO^^Mq have the cooling timescale comparable with the 2-body heating 
timescale. Thus, the relatively poor clusters w ith M < W^^Mq in Lbox = ISO/i-^Mpc models may 
suffer from the artificial 2-body heating. Since it is difficult to predict the possible systematic effect 
based on the timescale argument alone, it is most straightforward to quantify the effect by a careful 
comparison with the results in Lbox = 75/i~^Mpc runs which are almost free from the artificial 
2-body heating. In the following analysis, we use clusters with M > lO^^M© and M > 10^^ Mq for 
L150 and L75 models, respectively. Table 2 indicates the number of those clusters which satisfy 
the above criteria for each model at different redshifts. 

It should also be noted that one should not worry the numerical 2-body heating at the central 
regions of clusters where tdyn S> tcooi even if tcooi ^ hhody As noted in Steinmetz & White (1997), 
those regions will experience a catastrophic cooling; in fact, they are exactly the place where we 
attempt to suppress the overcooling by decoupling the cold gas particles. 
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3.2. Temperature — mass relation 

A conventional analytical modeling of clusters of galaxies assumes that the ICM is isothermal 
and in hydrostatic equilibrium within the dark matter potential. Then the ICM temperature is 
predicted to be 

kbTx = l—;r 

orvir 

. 5.2. f^!!^) \ U + '.d keV (7) 

in terms of the cluster mass M, where Ac is the mean density of a virialized object at a formation 
redshift Zf, and 7 is a fudge factor of an order unity; if the cluster is isothermal and its one- 
dimensional velocity dispersions is equal to \/ GM/3rvir, then 7 is the inverse of the spectroscopic 
/^-parameter and approximately given by 1.2 (Kitayama & Suto 1997). 

Since the above Tx-M relation is the most important ingredient in translating the Press - 
Schechter mass function into the XTF, the reliability of the conventional cluster abundance crucially 
depends on the applicability of this relation. Figure 5 plots the Tx-M relations for non-radiative 
and multiphase models at z = 2, 1 and 0. In order to avoid the possible artificial two-body heating 
effect, we show the results for clusters of M > 10^^ M© and M > 10^^ M© in L75 and L150 models, 
respectively. 

This figure implies three major conclusions; first, comparison of the upper and middle panels 
indicates that the current numerical resolution is sufficiently good and the Tx-M relation seems 
to be converged. Second, the mass weighted temperature is well fitted to equation (7) with 
7 = 1.2, while the emission weighted temperature Tx is systematically higher for the same M. 
Finally, the simulated Tx-M relation is almost unaffected by the phenomenological (in the current 
simulation) treatment of the thermal evolution of the ICM gas, and the results of L150A and 
L150UJ are almost identical. We also verify that the results of L150CJ axe almost identical to 
those of LI SOU J implying that the presence of a UV background does not affect our conclusions. 
This is the first successful check of the relation, made possible by the sufficiently large volume and 
good resolution of our simulations. 



3.3. Luminosity — temperature relation 

In contrast to the Tx-M relation, the luminosity - temperature relation of X-ray clusters is 
not easy to predict. This is because the X-ray luminosity is proportional to the square of the gas 
density, and thus sensitive to the thermal evolution of the ICM. A simple self-similar model, which 
predicts Lx (x rJ-(l-l-Zf)^/^ (Kaiser 1986), is shown to be completely inconsistent with the observed 
relation of Lx oc (1 + z)'^ where 2.6 < a < 3.0 and C ^ (Edge & Stewart 1991; David et al. 
1993; Markevitch 1998). 
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The Lx - Tx relation of our simulated clusters (Fig. 6) also exhibits the difficulty to obtain 
a reliable estimate of the luminosities. While it is reasonable that the results are sensitive to the 
ICM thermal evolution, they are also affected by the numerical resolution. The cluster luminosities 
in Lbox = 150/i~^Mpc models are systematically underestimated relative to those in 75/i~^Mpc 
models. Since we employ the equal-mass particle in the simulations, the resolution problem would 
be more serious for smaller clusters. Actually this explains why L150A model produces Lx oc T^- 
accidcntally although this non-radiative model should result in Lx oc r|-. Thus we conclude 
that the reliable estimate of the X-ray luminosities requires much higher resolution, especially for 
less luminous clusters, than ours in addition to the more physical implementation of the galaxy 
formation process. Therefore, we disagree with Pearce et al. (1999b) who concluded that the effect 
of cooling suppresses the X-ray luminosities; while our results do not reach convergence either, we 
find the opposite trend, i.e., the luminosities increase with the cooling (Figs. 1 and 6). After all, 
since their resolution is similar to ours, their results cannot escape from the resolution problem, 
and it is premature to draw any conclusion on the luminosities from cosmological SPH simulations. 



3.4. X-ray Temperature Function 

Finally we show the XTFs in non-radiative and multiphase models at rcdshift z = 0.0 and 1.0 
(Fig. 7) as a function of the emission and mass weighted temperatures. These should be compared 
with the analytical prediction on the basis of the Press-Schechter mass function (Press &: Schechter 
1974) and the Tx-M relation mentioned above. Since we have already showed that the Tx-M 
relation agrees well with the analytical expectation, our simulated clusters in cosmological volumes 
can examine the statistical reliability of the analytical prediction of XTF for the first time. Figure 7 
implies that the analytical and numerical results agree well with each other almost independently 
of the ICM thermal evolution model; if we adopt 7 ~ 1.6 and 1.2 for XTFs in terms of emission- 
and mass-weighted temperatures, respectively. It should be also noted that since the XTFs from 
different numerical resolutions are nearly the same within the statistical errors, we can state that 
the XTFs from our simulations do not suffer from any serious numerical artifacts. These results 
justify the use of a simple analytical model for the cluster abundance extensively applied before 
(White, Efstathiou & Frenk 1993; Viana & Liddle 1996; Eke, Cole & Frenk 1996; Kitayama & Suto 
1997; Kitayama, Sasaki k Suto 1998). The observed XTF of local {z < 0.1) clusters by Markevitch 
(1998) is also shown in Figure 7 and lower by factor of 2-3 than the simulated ones at Tx — 3-9 
[keV]. 



4. Conclusions 

On the basis of a series of the large cosmological SPH simulations, we have specifically addressed 
the reliability of the analytical predictions on the statistical properties of X-ray clusters. In order 
to distinguish numerical artifacts from real physics, we performed simulations with two different 
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numerical resolutions. Our main findings are summarized as follows: 

(i) The inclusion of radiative cooling in the high-resolution simulations substantially change 
the luminosity of simulated clusters. Without implementing the galaxy formation scheme, this 

leads to an artificial overcooling catastrophe as demonstrated by Suginohara & Ostriker (1998). 
With a phcnomenological prescription of cold gas decoupling like Pearce et al. (1999a), however, 
the overcooling is largely suppressed. Nevertheless the predicted X-ray luminosities of clusters are 
not reliable even by one order of magnitude. 

(ii) In contrast to the huge uncertainties on the X-ray luminosities, the temperature of simu- 
lated clusters is fairly robust both to the ICM thermal evolution and to the numerical resolution, 
and in fact are in good agreement with the analytic predictions commonly adopted. We also showed 
that the emission-weighted temperature is 1.3 times higher than the mass-weighted one. 

(iii) The analytical predictions for the X-ray temperature function translated from the Press 
- Schechter mass function are fairly accurate and reproduce the simulation results provided that 

the fudge factor 7 = 1.2(1.6) is adopted in the mass-weighted (emission-weighted) temperature - 
mass relation (7). The XTFs simulated in the cosmology adopted in this paper are slightly higher 
than the observed one of local clusters by Markcvitch (1998). 

Our simulations have made statistically unbiased synthetic catalogues of clusters of galaxies in 
virtue of their large simulation volume and sufficient numerical resolution. Thus, they can be the 
theoretical references which should be compared with the results from future cluster observations 
with X-ray satellites including XMM and Astro-E, and can contribute to probing cosmological 
parameters from observed cluster abundances in due course. Additional radiative processes such as 
heat conduction of ICM and energy feedback from supcrnovae explosions, which may considerably 
affect the physical properties of ICM, must be considered as future works in order to solve the 
discrepancy of L-T relation. In addition, since we consider only one fairly specific cosmological 
model, it is necessary to perform simulations under other cosmological models. 

We thank an anonymous referee for pointing out the importance of the two-body heating in the 
simulations. K.Y. and Y.P.J, gratefully acknowledge the fellowship from the Japan Society for the 
Promotion of Science. Numerical computations were carried out on VPP300/16R and VX/4R at 

the Astronomical Data Analysis Center (ADAC) of the National Astronomical Observatory, Japan, 
as well as at RESCEU (Research Center for the Early Universe, University of Tokyo) and KEK 
(High Energy Accelerator Research Organization, Japan). This research is supported by Grants- in- 
Aid by the Ministry of Education, Science, Sports and Culture of Japan to RESCEU (07CE2002), 
and by the Supercomputer Project (No. 99-52) of KEK. 
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Table 1. Summary of the simulation models 



model 


Lbox[Mpc] 


cooling 






mdark[M0]^ 


cold gas decoupling 


L150A 


150/i-i 


off 


0.0 


2.0 X IQi" 


1.7 X 10" 


No 


L150C 


150/1"^ 


on 


0.0 


2.0 X 10i° 


1.7 X lOii 


No 


L150UJ 




on 


1.0 


2.0 X IQio 


1.7 X 10" 


Yes 


L150CJ 




on 


0.0 


2.0 X IQio 


1.7 X 10" 


Yes 


L75A 


75/1-1 


off 


0.0 


2.4 X 10^ 


2.2 X 101° 


No 


L75C 


75/1-1 


on 


0.0 


2.4 X 10^ 


2.2 X 10i° 


No 


L75UJ 


75/1-1 


on 


1.0 


2.4 X 10^ 


2.2 X 10i° 


Yes 



gas and dark matter mass per particle 



Table 2. Number of simulated clusters that we used in the analysis 



model 


z = 2.0 


z = 1.0 


z = 0.0 


L150A 


12 


72 


224 


L150UJ 


11 


76 


220 


L150C 


12 


75 


222 


L75A 


123 


263 


295 


L75UJ 


117 


241 


272 


L75C 


120 


245 
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Fig. 1. — Spherically averaged profiles of dark matter and gas densities, ICM temperature and 
the integrated X-ray luminosity from the center, for a cluster at z = 0.0 (solid lines) and z = 1.0 
(dashed lines) in models L150A, L150UJ and L150C. Vertical lines indicate our lower limit on the 
SPH smoothing length, /imin- The viral mass of each cluster at z = 0.0 is quoted in the middle 
panels. For reference, the universal density profile at z = (Navarro, Frenk, & White 1997) 
corresponding to the virial mass of the cluster is plotted in open squares (top panels). 
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Fig. 2. — Same as Figure 1 but for models L75A, L75UJ, and L75C at z = 0.5 (solid lines) and 
z = 1.0 (dashed lines). The virial mass of each cluster at z = 0.5 is quoted in the middle panels. 
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Fig. 3. — Radial profiles of the local dynamical timescale, t^yn, the artificial two-body heating 
timescale t2body a^id the cooling timescale tcooi for three different clusters at 2: = 0.0 (lower panels) 
and z = 1.0 (upper panels) in model L150UJ. Their virial mass and temperature are indicated in 
each panel. 
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Fig. 4.— 



Same as Figure 3 but for model L75UJ. 
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Fig. 5. — Temperature-mass relations for non-radiative and multiphase models at z = 0.0, 1.0 and 
2.0. The lowest panels show the emission-weighted temperature while others the mass-weighted 
temperature. Lines indicate the theoretical model (eq.[7]) with 7 = 1.6 for emission- weighted 
temperature, and 7 = 1.2 for mass-weighted one. 
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Fig. 6. — Luminosity-temperature relations for non-radiative and multiphase models at z = 0.0, 
1.0 and 2.0. Thin lines indicate the relations predicted from the self-similar evolution. Bold lines 
indicate the empirical relation from observations (David et al. 1993). 



- 18 - 




Fig. 7. — X-ray temperature functions at z = 0.0 (lower panels) and z = 1.0 (upper panels) for 
each model. The left panels uses emission weighted temperature, whereas the right ones mass 
weighted temperature. Lines are theoretical predictions using Press-Schechter mass function and 
T-M relation (eq.[7]) with 7 = 1.0, 1.2 and 1.6. Shaded regions indicate the observed XTF of local 
{z < 0.1) clusters by Markevitch (1998). 



